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The Tavis-Cummings model for more than one qubit interacting with a common oscillator mode 
is extended beyond the rotating wave approximation (RWA). We explore the parameter regime in 
which the frequencies of the qubits are much smaller than the oscillator frequency and the coupling 
strength is allowed to be ultra-strong. The application of the adiabatic approximation, introduced 
by Irish, et al. (Phys. Rev. B 72, 195410 (2005)), for a single qubit system is extended to the 
multi-qubit case. For a two-qubit system, we identify three-state manifolds of close-lying dressed 
energy levels and obtain results for the dynamics of intra-manifold transitions that are incompatible 
with results from the familiar regime of the RWA. We exhibit features of two-qubit dynamics that 
are different from the single qubit case, including calculations of qubit-qubit entanglement. Both 
number state and coherent state preparations are considered, and we derive analytical formulas that 
simplify the interpretation of numerical calculations. Expressions for individual collapse and revival 
signals of both population and entanglement are derived. 
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I. INTRODUCTION 

Two level systems that interact with a harmonic oscil- 
lator can model many physical phenomena, such as nu- 
clear spins interacting with magnetic field atoms in- 
teracting with electromagnetic field 0, [|| , electrons cou- 
pled to a phonon mode of a crystal lattice super- 
conducting qubits interacting with a nano-mechanical 
resonator a transmission line resonator @, Q, or 

an LC circuit [3, etc. The dynamics of all such sys- 
tems is governed by the Rabi Hamiltonian [l[ : 

jjRabi = + ^ a t a + hu^(a + a)){a+ + Sr_), (1) 

where the a z and a + + <r_ — a x are the usual Pauli ma- 
trices in the Hilbert space of the qubit and &' and a refer 
to the creation and annihilation operators of an inter- 
acting mode of a harmonic oscillator. Although studied 
extensively since it was first introduced in the context of 
nuclear magnetic spin resonance, analytical solutions for 
the eigenvalues and eigenfunctions of the Rabi Hamilto- 
nian still do not exist. 

In physical situations where the qubits are nearly res- 
onant with the oscillator and the coupling strengths be- 
tween the qubits and the oscillator are much smaller 
than the qubit and the oscillator frequencies, it is a 
good approximation to drop the counter rotating terms: 
a)(j + and a<r_, from ([T]) to obtain the so-called Jaynes- 
Cummings (JC) model with the Hamiltonian [2j: 

ftjc = .^a- z + ^a t a+faj^(a<7 + + CT_a + ). (2) 
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Under this approximation, called the rotating wave ap- 
proximation (RWA) , the dynamics of the system can be 
obtained in closed form @, 0] ■ 

A generalization of the JC model, called the Tavis- 
Cummings (TC) model, was introduced in the context of 
quantum optics to describe the collective behavior of mul- 
tiple atomic dipoles interacting with an electromagnetic 
field mode [lll - fl3j . The TC model has gained renewed 
interest as it can be used to implement quantum infor- 
mation protocols with the oscillator transferring informa- 
tion coherently between qubits [hi ]. Intrinsically multi- 
qubit properties such as quantum entanglement can be 
explored with the TC model in a variety of ways, employ- 
ing various entanglement measures such as concurrence 
for mixed-state pairs of qubits [l5j |. quantum negativity 
for slightly larger systems [l6j], and Schmidt weights for 
bipartitions of arbitrarily dimensioned pure multi-qubit 
states [I?] ]. 



II. MULTI-QUBIT BREAKDOWN OF THE RWA 

With recent advances in the area of circuit QED, it 
is now possible to engineer systems for which the qubits 
are coupled to the oscillator so strongly, or are so far 
detuned from the oscillator, that the RWA cannot be 
used to describe the system's evolution correctly [Tcj - Eol ] . 
The parameter regime for which the coupling strength is 
strong enough to invalidate the RWA is called the ultra- 
strong coupling re gim e (21rl27| . Niemczyk, et al. and 
Forn-Dfaz, et al. [19| have been able to experimentally 
achieve ultra-strong coupling strengths and have demon- 
strated the breakdown of the RWA. Motivated by these 
experimental developments and the importance of under- 
standing collective quantum behavior, we investigate a 
two-qubit TC model beyond the validity regime of RWA. 
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FIG. 1: Diagrams showing energy-level configurations: (a) 
compatible with the RWA, A <C u>o ~ w; (b) incompatible 
with the RWA, A > uo C The states \e) and |<?) are the 
eigenfunctions of a z : a z \e) = |e) and o" z |<?) = —\g). 



are the first analytic expressions for the individual revival 
signals beyond the RWA, as well as analytic expression 
for the collapse and revival dynamics of entanglement. In 
the quasi-degenerate regime, the invalidity of the RWA in 
predicting the dynamical evolution will clearly be demon- 
strated in Sec. |V](see Figs. IHandE]). 

We begin with a generalization of ([I} in which the a 
operators are replaced by two-qubit counterparts [l^ : 



H = hui^Sz + Hud/ a + hujf3{a + a )S X 



(3) 



where 



S z = + ai 2 >), and S x = + a^). (4) 

In experiments dealing with artificial qubits, such as 
Cooper pair boxes, it is possible to bias the qubits, which 
results in an additional term in the Hamiltonian: 



The regime of parameters we will be concerned with is the 
regime where the qubits are quasi-degenerate, i.e., with 
frequencies much smaller than the oscillator frequency, 
luq <C uj, while the coupling between the qubits and the 
oscillator is allowed to be an appreciable fraction of the 
oscillator frequency. In this parameter regime, the dy- 
namics of the system can neither be correctly described 
under the RWA, nor can the effects of the counter rotat- 
ing terms be taken as a perturbative correction to the 
dynamics predicted within the RWA by including higher 
powers of j3. For illustration, systems are shown in Fig. 
[T] for which the RWA is valid, or breaks down, because 
the condition ujq w uj is valid, or is violated. The regime 
that we will be interested in, for which loq <C uj, is shown 
on the right. 

Prior numerical work by Irish has been directed to the 
dynamics of a single quasi-degenerate qubit interacting 
with an oscillator in the ultra-strong coupling regime, 
and carried out by developing an adiabatic approxima- 
tion [2l|, with an extension to a generalized RWA [22j], 
and also by Hausinger, et al., by using van Vleck pertur- 
bation theory (26j . The adiabatic approximation and van 
Vleck perturbation theory were shown to work best for 
small qubit frequencies and high coupling strengths. The 
adiabatic approximation was shown to fail in the regime 
where the JC model works well, i.e., when the qubit is 
resonant with the oscillator and the coupling is small. 
This gap between the regime of validity of the adiabatic 
approximation and the regime of validity of JC model 
was bridged by the generalized RWA 22j, which works 
well in both regimes. 

Here, within the adiabatic approximation, we extend 
the examination to the two-qubit case. Qualitative differ- 
ences between the single-qubit and the multi-qubit cases 
are highlighted. In particular, we study the collapse and 
revival of joint properties of both the qubits. Entan- 
glement properties of the system are investigated and it 
is shown that the entanglement between the qubits also 
exhibits collapse and revival. We derive what we believe 



Hbias — heS x , 



(5) 



where e is called the static bias. Taking finite bias into 
account, e ^ 0, an analysis of the dynamics of a single 
qubit interacting with a harmonic oscillator beyond the 
RWA in the ultra-strong coupling regime was done in 
25, 261. Here we assume that e = 0. 



III. INFORMAL ANALYSIS 

Before proceeding with a detailed treatment, we note 
that an informal approach to the Hamiltonian ([3]) is pos- 
sible, and can be helpful in interpreting further analysis. 
The disparity in time scales signaled by the inequality 
ujq <C uj governs new effects that will occur. To see this, 
we let uio be sufficiently small as to be negligible, thus 
removing the S z from any role in H. Then S x becomes 
constant, say S x (0). The Heisenberg equation for the re- 
sponse of the oscillator amplitude a becomes trivial, with 
the solution 



a(t)+l3S x (0) = (a(0)+j8£ x (0))e- 



(6) 



which is easily interpreted in the expected-value sense: 
the evolution of (a) is sinusoidal at frequency w and cen- 
tered at -P(S x (0)). 

Of course, S x is not constant if ujq ^ 0. As an operator, 
it evolves in time. Its evolution is determined by the 
commutator with H , and this leads to the three coupled 
Bloch-type equations: 



dS x /dt = -u) S y , 
dSy/dt — uj S x — /3ui(a - 
dS z /dt — [iuj(a + aJ)S y 



(7) 

(8) 
(9) 



Under realistic current laboratory conditions (3 -C I, so 
unless the oscillator amplitude is very great the main spin 
motion is a slow precession of S x and S y at frequency 
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wo, with small and very rapid oscillations at frequency oj 
arising from the f3(a + a^) terms. 

Two comments are obvious at this level of analysis. 
First, since S x changes nearly periodically on the time 
scale ~ 27t/ojo, we expect the center of oscillator motion 
to follow these slow changes back and forth. Second, 
the rapid oscillations around the slow precession are of 
both signs ±w, so they contain the effect of the counter- 
rotating terms omitted by the JC model. It should be 
noted that this informal analysis is not specific to any 
particular number of qubits and all the comments of this 
section are equally applicable to an if-qubit system. 



IV. SPECTRUM OF H 

We first find the eigenspectrum of H when ujq = 0. 
The Hamiltonian without the RWA then takes the form 



Hq = fvjja/a + h/3tj(a + a/)S x 



(10) 



The eigenstates and eigenvalues of Hq satisfy the eigen- 
value equation: 



a)a + P{a + a f )S x |$) = E\<S> 



(11) 



The eigenstates | $) will be products of qubits and oscil- 
lator states, and take the form 



I*) = \ j,m)\4>„ 



(12) 



Here \j, m) are the eigenstates of S x and \4> m ) are the 
oscillator eigenstates found from Hq by replacing S x by 
its eigenvalue corresponding to \j, m) [2l| . 
The four eigenstates of 



\j,m) = |1,±1), |1,0) and |0,0), 



(13) 



with eigenvalues m. In terms of the simultaneous eigen- 
states of ai and dx , Sx |±) = ±|±), the states \j,m) 
can be written as: 




/l 0\ 

o 1/V2 1/V2 

1/V2 -1/V2 

\0 1/ 




(14) 



Having found \j, m), let us now find \4>m) that satisfy 
the eigenvalue equation: 



hjj [at a + mj3(a + at)] \<j) m ) = E\<j) n 



(15) 



We denote m/3 by /3 m , which we take real. Then by com- 
pleting the square in (|15p , we get a new number operator 
equation: 

(a t +/3 m )(a + /3 m )|</> m ) = (E/tuv + 0^) \<f> m ) 

= N\4> m ), N = 0,1,.... (16) 



Using the displacement operator, D(a) = exp[a(a^ — a)] 
(for real a), we can write the expression on the left side 
of (|TH)) as (p m )at aD(p m )\<Pm) ■ Then multiplication of 
this by D(j3 m ) converts (fTBf into 



a) a 



(D{p m )\4> m )) =N(D(j3 m )\<p m )), (17) 



which shows that the original oscillator and its displaced 
counterpart have the same eigenvalues, and relates their 
eigenstates as 

D(fi m )\<t%) = \N) or 

|0£) = D{-p m )\N) = \N m ). (18) 

Thus, finally, the joint qubit-oscillator eigenstates are of 
the form: 

I*) -> \$j,m,N) = \j,m) \N m ), (19) 
and the energy E in ([TB]) takes values: 

E Ntm = fiu(N-^). (20) 

Thus, we see that depending upon the state of the 
qubits, determined by \j,m), we have four harmonic 
oscillator potential wells in x — p phase space, where 
x = [a) + a)/y/2 and p = i{a) — a)/y/2. These poten- 
tial wells have their equilibrium positions displaced by an 
amount proportional to 2m/3. For m — 0, the oscillator 
potentials are not displaced, whereas for m — ±1, they 
are displaced in equal and opposite directions. A very 
important thing to note from (|2U| is that the eigenstates 
with the same value of N are not degenerate, e.g., the 
states |l,0)|Ao) and |l,l)|JVi) differ in energy by fuv/3 2 . 
For contrast, in the single qubit case, when wo = 0, the 
eigenstates of the Hamiltonian with the same value of N 
remain degenerate irrespective of the value of (3 21( . 

The three potential wells corresponding to the states 
\l,m)\N m ) are schematically shown in Fig. [2j The dis- 
placement of the equilibrium position of the potential 
wells and the relative lowering of the energy levels for 
m = ±1 states is evident from the figure. 

One may say that because of its coupling to the qubits 
the original oscillator is not really the "effective" oscilla- 
tor, with the consequence that a definite number of its 
excitations does not correspond to a definite number of 
the effective excitations, and vice versa. This is the na- 
ture of the displacement operation. In a discussion of 
two level systems interacting with a harmonic oscillator 
beyond the RWA, the use of a displaced harmonic oscilla- 
tor basis was first used by Schweber 28]. The displaced 
oscillator states have the properties: 



(N m \N^) = 5 NtN ,, 
(N m \N^) ? 0, 

and in particular 

(Nx\N ) =e-^' 2 L N {p 2 ), 



(21) 



(22) 
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m=0 

m=l m=-l 




|l,±l)|7V±i), are nearly degenerate. We call this quasi- 
degenerate triplet of states the N th manifold. Because of 
finite u>q, there will be transitions between various states: 
|l,m)|AT m ) and \l,m')\N m ,). These transitions can be 
classified under two categories: (a) transitions that take 
place between levels belonging to different manifolds and 
(b) transitions that take place between the three states 
that belong to the same manifold. For two adjacent man- 
ifolds, transitions of type (a) are shown in Fig. [3][a) and 
for the same two manifolds, transitions of type (b) are 
shown in Fig. 02[b). Suppose the states \l,m}\N m ) and 



FIG. 2: The three potential wells corresponding to the states 
|l,l)LVi) (left), |1,0)|JV > (middle) and |l,-l)LV_i> (right). 
The factor AX zp is the zero point fluctuation of a harmonic 
oscillator. For an oscillator of mass M and frequency ui the 
zero point fluctuation is given by AX zp = ^/h/2Muj. 



where Ln(x) is a Laguerre polynomial. The non- 
orthogonality condition, (N m \N' m ,) 7^ 0, plays an impor- 
tant role in subsequent analysis. 

Next, we extend the discussion to examine the eigen- 
spectrum of H when uj n ^ 0. Using the basis |j, m)\N m ) 
we now look for the eigenstates and eigenvalues of H 
when w 7^ 0. We note that because |0, 0) is a simultane- 
ous eigenstate of S z and S x : 



(a) 



&|0,0) = 0, 

s x \o,o) = 0, 



(23) 



the states |0, 0) | AT ) (for any N) are eigenstates of H 
with eigenvalues £jv,o = fiNu, even when luq 7^ 0. This 
allows one to find the exact evolution of the system in the 
projected Hilbert space spanned by the states |0, 0)|iVo). 
However, finding the evolution of a state spanned by 
|l,m)|AT m ) is a challenge because the states \l,m)\N m ) 
are not simultaneous eigenstates of S z and S x . We now 
look at the Hamiltonian that is spanned by \l,m)\N m ) 
states only. We assume that even though wo 7^ 0, it is 
still small compared to the frequency of the oscillator 
as a result of which one can treat the term HujqSz as a 
perturbation to the energy spectrum found for the case 
when luq = 0. In particular, we will restrict our anal- 
ysis to the regime: uio < 0.25w, which we label as the 
quasi-degenerate regime. 

We start by noticing from (f20| that when uj n = and 
/3 2 is close to an integer, say p, the three states: 1 1, 0) | N ) 
and |1, ±1)| (AT + p)±i), are grouped together in energy 
and are nearly degenerate. In what follows, we will not 
be concerned with very high values of |/3|, but explore 
the regime that is experimentally achievable with cur- 
rent technology or is likely to be realizable within the 
near future. For this reason, we restrict our analysis to 
the regime where \(3\ < 0.25, which is strong enough to 
invalidate the RWA, i.e., which lies in the ultra-strong 
coupling regime. Under this assumption, states with 
the same value of oscillator excitation: |l,0)|Ao) and 
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FIG. 3: (a) Transitions induced by hw^Sz between states of 
different manifolds, (b) Transitions induced by huoS z be- 
tween states of the same manifold. 



|1, m')|iV m /) belong to different manifolds, 
tion matrix element between them is 



w (l,m\S z \l,m')(N m \N' n 



The transi- 



(24) 



If the above transition matrix element is much smaller 
than the energy difference between them, i.e. 



w 



[l,m\S z \l,m'){N m \N m ,) 4Zu\N-N'\, (25) 



then the transitions of type (a) would be energetically 
suppressed. On the other hand, because transitions of 
type (b) occur between nearly degenerate states, there 
could be appreciable transfer of population between 
them. Based on these arguments, one can neglect all ma- 
trix elements in H that lead to transitions between dif- 
ferent manifolds and retain only those terms that induce 
transitions between states of the same manifold [U [29l— 
l3l| . This approximation was used by Irish et al. |2l| 
to study the dynamics of a single quasi-degenerate qubit 
interacting with a high frequency oscillator. 
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Under the above assumption of neglecting transitions 
between states belonging to different manifolds, the 
Hamiltonian spanned by the |1, m)\N m ) basis reduces to 
3x3 block diagonal form with each block corresponding 
to a given manifold. For the N th manifold, the Hamilto- 
nian takes the form 



N — Q N 

n N n n N 

n N N-/3 2 



(26) 



where rows and columns are arranged in the order: 
|1, l>]7Vi>, |1,Q)|JV ) and |1, -l)|iV_i), and the off- 
diagonal terms are the normalized Rabi frequencies given 
by: 



N = ^(1,11^11, Q)(JVi|JV ) 



2 uj 



(27) 



In writing Hn, we have used the fact that (iVi|JVo) is real 
and is equal to (iV_i|iVo). Note that for any value of N, 

We note from the form of Hn that due to the pres- 
ence of the off diagonal elements, the state of the oscilla- 
tor, which is displaced depending upon the state of the 
qubits, changes with the changing state of the qubits. 
This change happens in a time scale of 1/ (luD,n) > 1/luq 
which is much slower than the characteristic time scale 
of the oscillator, which is 1/lu. One thus sees that the 
oscillator state adiabatically adjusts itself to the state of 
the qubits. For this reason, the abov e appro ximation is 
known as "adiabatic approximation" [2l|, [31| • 

The unnormalized eigenfunctions and eigenvalues of 
Hn are: 



\£°n) 




'N 



V i 

hw{N-p 2 ), 



V. POPULATION DYNAMICS 

Analysis of the dynamical properties of a single qubit 
in the RWA-violating quasi-degenerate regime can be 
found in [2l|, [32| . Here we take a step in the direction 
of A'-qubit evolution by considering A=2, and defer an 



introduction to cases for K > 2 to Sec. I VIII We will 
stay within the quasi-degenerate parameter regime men- 
tioned in Sec. IIVI (cjq < 0.25w) and further make the 
following assumptions that would considerably simplify 
the expressions in (f28|) : 



n N >/3 2 



< 0.2 and N > 0. 



(29) 



This allows some obvious simplifications: 8il N — /3 4 
8il N and N — (3 2 w N respectively. Then the expressions 
for the eigenfunctions and eigenvalues of Hn simplify to: 



14) 




tvuiN, 

hco (n± V2fijv) • 



(30) 



For illustration, consider initial states that belongs to 
the N th manifold: 



|*±(0)) = |1,±1)|7V ±1 ), 



N, 



± T2 K) 



(31) 



Using ([3^]). the probability amplitude for the qubit to 
remain in the initial state is easily found to be 



(*±(0)|*±(t)) 



I + cos(V2~n N ujt)). (32) 



When squared, the probability shows two frequencies 
of oscillation, V^jvw and 2v / 2^ato;. Since three basis 
states are involved, we could expect three frequencies, 



N 



This is in 



but two are equal: 
contrast to the single-qubit case where only one Rabi 
frequency determines the evolution [2l|, |32[ . The two fre- 
quencies contribute to the probability as follows: 

3 1 1 

P lt±1 (N,t) = - + - cos(V2Q N ujt) + - coa(2V2n N u}t). 
8 2 8 

(33) 

A different initial state, also characteristic of the 
two-qubit case, that belongs to the N th manifold is 
(28) |1,0)|A ) = (\£%) + \£ N ))/V2. One finds that the prob- 
ability to remain in this state oscillates with only one 
frequency, but twice as great as the higher frequency in 
the previous example: 



P h0 {N, t) = l - + 1 cos(4V2~n N cjt). (34) 

A number state is in most cases not a reasonable model 
for describing experimental excitation of the oscillator. A 
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coherent-state description for the oscillator is more real- 
istic, and in that case the Poisson distribution of num- 
ber states creates a distribution of Rabi frequencies Qnoj. 
The oscillations of solutions for different N values rapidly 
get out of phase with each other and the signal collapses 
quickly. However, the main contribution of the Poisson 
distribution comes from its peak near N = n w \a\ 2 
where adjacent Rabi frequencies differ by a small com- 
mon amount 8£l(n)u), which leads to a rephasing of the 
main terms in the summation at integer multiples of the 
time 2tt / Sfl(n)uj . One thus expects a sequence of revivals 
and then re-collapses of the signal, which are familiar in 
parameter regimes where the RWA is valid (33j . 

The collapse and revival behavior in the adiabatic ap- 
proximation for a single qubit case was studied by Irish 
et al. [2l[ and Sandu [32[ and we explore here the two- 
qubit counterpart. We obtain analytical expressions for 
the collapse and revival times and also for the individual 
revival signals. Since in the two-qubit case the system 
eigenstates are displaced number states, a coherent-state 
sum of them produces a displacement of the coherent 
state | a) as the initial state: 



|*_(0)) = |l»-l)l«-i> 

= \1,-1) D(/3)\a). 



(35) 



Then the probability for the qubits to remain in the state 
|1, — 1) is found to be: 

Pi,-x(a,t) = ^ + ±S(t,uj ) + ~S(t,2LJ ), (36) 



where 



S(t,u ) = i^r 1 — coB(w <JVi|JVo)i). (37) 



JV=0 



Nl 



If the average excitation of the oscillator, \a\ 2 , is large, 
one can evaluate the above sum approximately (see Ap- 
pendix) to get 



S(t, wq) = Re 



.fc=0 



(38) 



where 



Re [Sk(t,u>o)] = exp 



2(l + (7Tfc/) 2 ) 
COS ($ im ) 



(1 + (TTkfY) 



2^/4- 



(39) 



In (f3"9"| we have defined 



r = ^te-^'\ 
f = l«/?| 2 , 

r fe = 2^fc(l + //2)//3 2 , 



(40) 



and $/ m is given in (|A.10|) . 



From (|38|) and (|39|) . it is clear that S(t,uo) exhibits 
collapse and revival with Sk{t,u>o) describing the evolu- 
tion around the k th revival time. These individual revival 
signals, Sk(t, too), have three salient features: (a) the ex- 
ponential term in (J39J) determines the envelope of the 
revival signal, (b) the cosine term governs the fast os- 
cillatory dynamics and (c) the factor in the denominator 
determines the height of the k th revival. The revival time 
and the height of the k th revival are: 



2nk 



/;■"■ = — -(l + |a/3| 2 /2), 
1 



k ~^/3 2 ' 



hu = 



(l + fc 2 ^ 2 |a/3| 4 ) 1/4 ' 



(41) 



As usual, the revivals are periodic and the heights of the 
revivals successively decrease and thus the revivals are 
never complete. The revival time increases with the in- 
crease of the oscillator excitation amplitude, \a\, and de- 
creases if the coupling parameter, \/3\, is increased. From 
(1391) . we note that the width of the primary revival, for 
which k = 0, is 



5t = 



1 



and the width of the k th revival is given by 



5n = 5T ^l + (7rk\ap\ 2 y 



(42) 



(43) 



Thus, we see that the width of the successive revival sig- 
nals keep increasing. We note that the term \a(3\ 2 /2 in 
the expression for the revival time (|4ip. is an improve- 
ment over the results given by Irish et al. [2l| and Sandu 

From (I36|) we note that two functions are responsible 
for the evolution of Pi ) _i(a,t): S(t,uj ) and S(t,2uJo). 
Thus, we get two different revival sequences in the evo- 
lution of Pi ) _i(a, t). The analytic formula derived for 
Pi,-i{ct,t) is plotted in Fig. [4] and is compared with the 
numerical calculations. The revival signals correspond- 
ing to the terms Sk(t,u)o) and S2k{t, 2ujq) overlap in time 
and produce a beat note. This is evident in Fig. @] As 
mentioned earlier, the RWA completely breaks down in 
the parameter regime we consider. This can clearly be 
seen in Fig. 0] where the evolution of Pi ; _i(o;, t) calcu- 
lated within the RWA disagrees even qualitatively with 
the numerical calculations. 

The expression (j39|) , was derived under the constraint 
| a/3 1 <C 1 (see (|A.3I) 1 . Within this constraint, the revival 
time was found to be a monotonically increasing and de- 
creasing function of \a\ and \/3\ respectively. If the oscilla- 
tor excitation number and the coupling strength are not 
restricted by the constraint \a/3\ -C 1, the revival time is 
no longer a monotonic function of \a\ and This non- 
monotonic behavior was numerically explored in [2l[ . In 
the limit of very high oscillator excitation number, we 
can employ the asymptotic expression for L„(/3 2 ) to de- 
rive an analytic expression for the revival time of S(t, ujq) 
that is not restricted by the constraint \a/3\ <C 1. 
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FIG. 4: Collapse and revival dynamics for Pi j _i(q, t), given 
ujo = 0.15a;, /3 = 0.16 and a — 3. Note the breakup in 
the main revival peaks of the numerical evaluation, which 
comes from the loq-2ujo beat note, not included in the analytic 
calculation. The RWA is seen to break down completely in 
the parameter regime considered. 



As we already mentioned, revivals should occur at mul- 
tiples of the time t = t rev such that 

5Cl(n)u}f ev = 2tt, or 

Wo e- /32/2 |^+i(/3 2 ) - M/3 2 )| t rev = 2tt. (44) 

If the oscillator is highly excited, n 3> 1, one can use 
the following asymptotic formula for the Laguerre poly- 
nomial 1341: 



Jim e- x / 2 L n (x) = 



cos (2\/nx — 7r/4) 
Y / 7r(nx) 1 / 4 



(45) 



to obtain: 

'u: t rev 
2tt 



cos (2\a(3\ - tt/4) 



+ 



J^l sin (2|a^|- tt/4) 



(46) 



From (|46p. the non-monotonic dependence on a and /3 of 
the revival time is clear. Note that equation f4"6"]) is not 
restricted by the constraint |a/3| -C 1. 

In Fig. [5J we plot S(t, ujq) for a = 10 and various values 
of (3. The revival times predicted by (|4"oT) are denoted by 
vertical lines and are seen to have excellent agreement 
with the numerically evaluated revival signals despite 
a strongly varying location of the revivals. The figure 
clearly demonstrates the non-monotonic dependence of 
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FIG. 5: Numerical evaluation of S(t, Wo) for a = 10 and vari- 
ous values of /3. Note that the revival time is a non-monoitonic 
function of /3. The vertical lines correspond to the revival 
time predicted by Eq. (|46[) and are seen to coincide with the 
numerically evaluated revival signals. 



the revival time on the coupling strength and also high- 
lights the departure from the formula for t rev derived in 
Sec. IVl under the constraint \aj3\ <C 1. The revival en- 
velopes seen in Fig. [5] are not approximately gaussians as 
was the case for the revivals studied in Sec. |V] A detailed 
discussion of the non-trivial structure of the revivals for 
big values of |a| and |/3| can be found in (2lj . 

We now contrast the dynamical evolution of the two- 
qubit TC model with the single qubit system. We assume 
that the initial state of the single qubit system is: 



|*.(0)) = |l/2,-l/2)|a_ 1/3 ), 



(47) 



where a x \l/2, -1/2) = -|l/2,-l/2) and £>(/3/2)|a) = 
|q:_ 1 / 2 ) - The state |^ s (0)) is a single qubit counter- 
part of the two-qubit initial state, |\&_(0)) given in ([35]) . 
in the sense that both the states are such that (a) the 
qubit(s) and the oscillator are initially uncorrelated, (b) 
the qubit (s) state is an eigenstate of S x with the low- 
est possible eigenvalue of m and (c) the oscillator is in 
a displaced coherent state. Under the adiabatic approx- 
imation, the evolution of the state can be analytically 
derived [2l[ and one finds that the probability for the 
qubit to be in the state 1 1/2, —1/2) evolves as: 



I2JV 



Pi,_i (M) =1/2 1+5] e-H cos ^fijvt) 

V N=0 ' / 

= l/2(l + S(t, wo)). (48) 

There is only one revival sequence for the single qubit 
system as a consequence of having only one Rabi fre- 
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FIG. 6: Collapse and revival dynamics for Pi _i(a,t), given 

ujo = 0.15oj, /3 = 0.16 and a — 3. Note the single revival 
sequence. Also, note that there are no breakups in the revival 
peaks in contrast to the two-qubit case (Fig. The RWA 
fails to describe the dynamical evolution even for the single 
qubit case. 



cillator in an undisplaced coherent state \a): 

\m) = ^(i++)+i--»h, 

= + (49) 

We continue with the approximations used to study the 
evolution of Pi i _i(a, t) in the previous section. In partic- 
ular, we focus on the changes arising from the presence 
of four effective oscillators. Then, given N 3> j3, we make 
the following approximation: 



\N )xiD(TP)\N ) = \N ±1 ), 



(50) 



which leads to 



|C(0)) = - ? =(|l,l)|ai) + |l,-l)|a-i)). (51) 
Using (f5U|) . one evaluates the time evolved state to be: 



Ear 



2V2 ^y/M 



t e -iS%t/h _|_ e -i£ N t/h\ 



x (|l,l)|iVi) + |1,-1)|JV_!)) 
+V2(e~ t£ " t/h - e- l£ « t/h )\l,0)\N o ) 

and from (|50p one sees that this reduces to: 



quency in the single qubit case. The analytic and nu- 
merically exact evolution of Pi _i (a, t) is plotted in Fig. 
The single revival sequence is evident from the figure. 
A discussion on the multiple revival sequences for the K- 
qubit TC model, within the parameter regime where the 
RWA is valid, can be found in 13511. 



VI. ENTANGLEMENT DYNAMICS 



The evolution of entanglement between several non- 
interacting qubits coupled to a single mode or many os- 
cillator modes, which act like quantum buses mediating 
information between the qubits, has been studied exten- 
sively (e.g., see [36j - T38j ) . In all these cases, the interac- 
tion between the qubits and the oscillator mode(s) were 
treated within the RWA. New time scales and qualita- 
tively new features arise in regimes where the RWA is 
invalid, mandating an extension of these previous results 
[39U42I ] . Applications may be important in areas of quan- 
tum information processing, where coherent control of 
entanglement may be vital. 

A generic illustration can start in a configuration with- 
out correlation between the oscillator and qubits. We 
place the qubits in one of the a x Bell states and the os- 



!£(*)> 



-M 2 /2 



E 



2V2 ^ Q VNl 



( e -*£p/* + e -*fi**/A) 1) + |i, -1)) 
y/2( e - t£ P/ h - e - 4£ «*/ ri )|l,0)l |JV ). 



This is particularly compact in the a z eigenbasis (<3z|e) 
|e), cr z \g) = — Is)), where it becomes: 



m) = —t^ E 



a 



N 



V2 ^ Q VN\ 



x (e-'^^lee) + e- i£ » t/h \gg)} \N ). (52) 

In order to study the entanglement dynamics between 
the two qubits, we first trace out the oscillator degrees of 
freedom to get the reduced density matrix for the qubits: 



^\t) = j2(N a \m)m\No), 



Iss) (ssl 



l(\ee)(e 

1 & 

+ o(E 5fc ^' 2w o)l53)(ee|+H.c. 



(53) 



fc=0 



At time t — 0, the state of the qubits is pure, but as 
time evolves, the reduced state of the qubits becomes 
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mixed. One can use concurrence to quantify entangle- 
ment between two qubits that are in an arbitrary mixed 
state 15]. Concurrence varies in the range from zero 
to one with zero denoting no entanglement and one de- 
noting maximum entanglement between the qubits. The 
density matrix, p^ 1,2 ^(t), is an example of a so called X- 
matrix [43]. Calculating the concurrence for an X- matrix 
is particularly easy and for p^ 1,2 \i) it is evaluated to be: 



C(t) = 



(54) 



As discussed in Sec. |V] this expression has periodic re- 
vivals with each term, Sk{t, 2ujq), in the sum centered at 
the k th revival. When they are well resolved and don't 
overlap and one is only interested in the envelope of the 
revivals, one can neglect the interference between the var- 
ious terms in the sum and an approximate expression for 
the concurrence is found to be: 



C(t)«5^|5 fc (t,2wo)|, 



fc=0 



E 



(TTfc/) 2 ) 174 

-(2r-r fc ) 2 |q/3 2 | 2 N 
2(1 + W) 2 ) . 



x exp 



(55) 



This expression for concurrence is plotted in Fig. We 
see that the entanglement between the qubits exhibits 
collapse and revival and the analytic formula agrees well 
with the envelope of the numerically evaluated result, 
predicting correctly the time, height and width of the 
individual entanglement revival signals. 
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FIG. 7: Numerical and analytical evaluation of the entan- 
glement dynamics between the two qubits for cjo = 0.15a;, 
f3 = 0.16 and a = 3. Entanglement between the qubits ex- 
hibits collapse and revival. The analytic expression agrees 
well with the envelope of the numerically evaluated entangle- 
ment evolution. 



VII. GENERALIZATION TO A'-QUBIT SYSTEM 

The analysis presented so far is restricted to the two 
qubits case. In this section, we will qualitatively sketch 



the procedure for extending the formalism of studying 
two quasi-degenerate qubits interacting with a high fre- 
quency oscillator to the A"-qubit system. The Hamilto- 
nian governing the dynamics of the A-qubit TC model 
is an obvious generalization of ([1]), where a operators are 
replaced by their A-qubit counterparts: 



where 



1 K 

i=l 



1 K 

W, and ^ = -^<t 

i=i 



(56) 



(57) 



The eigenfunctions of H when uq = are \j, m)|$ m ), 
where the states \j, m) are eigenfunctions of S x and the 
states |3> m ) = \N m ) are the generalization of the dis- 
placed Fock states defined in (fT51) . For simplicity in no- 
tation we will take A to be even, and then we have 



j = 0, 1, . . . , A/2, with m=-j,..., j. 



(58) 



The eigenvalue of H corresponding to the state 
\j,m)\N m ) is the same as before, -EV.m = Hlj(N — 
m 2 /3 2 ), which is independent of A. If we assume that 
(A/2) 2 /? 2 <C 1, the states with the same value of oscil- 
lator excitation number, A, will have nearly the same 
energy and will form a quasi-degenerate manifold. For 
a given j, a quasi-degenerate manifold consists of 2j + 1 
states (one state for each to). 

The matrix elements of the perturbation term, ujqS z , 
in the \j,m)\N m ) basis are: 

S 

uo{jMSz\j',m'){N m \N , m ,) = u J f{N m \N^) 



A /(j-TO)(j+TO+l) 



i ^ m, ' — l,m 



(59) 



We see from the above formula that the perturbation 
Hamiltonian connects only the states with the same value 
of total spin j. Under the adiabatic approximation, we 
retain only those terms of the perturbation Hamiltonian 
for which N — N'. Thus the Hamiltonian for the A- 
qubit case breaks into (A/2) + 1 disjoint Hamiltonians 
with each disjoint Hamiltonian corresponding to a given 
j. This is a consequence of having 5jj> in equation (|59[) . 
Furthermore, for a given j, the Hamiltonian assumes a 
block diagonal form in the \j, m)\N m ) basis, with each 
block corresponding to a particular value of N . Each of 
these blocks has a dimension of (2j + 1) x (2j + 1). This 
block diagonalization of the Hamiltonian is a consequence 
of the adiabatic approximation. Finding the eigenvalues 
and eigenfunctions of each of these block diagonal ma- 
trices allows one to study the dynamical evolution of the 
system analytically. 

Entanglement dynamics of the A-qubit TC model is 
important in understanding multi-particle quantum co- 
herences. Multipartite entanglement is still largely mys- 
terious, in the sense that no approach is known that pro- 
vides both necessary and sufficient criteria for arbitrary 
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mixed-state entanglement of more than two parties or 
even for just two parties if their Hilbert states have di- 
mensions greater than 2x3 

The pure-state situation is clearer. In that case one can 
use the so called Schmidt weight to reliably quantify en- 
tanglement between two-party states of arbitrary dimen- 
sions [l7| . Thus, if our initial state is pure, we can con- 
ceptually partition the composite system into two parties 
and study the entanglement dynamics between them by 
using Schmidt analysis. Note that our composite system 
consists of K + 1 parties (K qubits and an oscillator) and 
so there can be 2 A — 1 different bi-partitions. 



VIII. VALIDITY REGION 

For better appreciation of the zones of validity of the 
different approximate approaches to the RWA and quasi- 
degenerate evolution dynamics, we show a 3D represen- 
tation in Fig. [8] The three axes in the figure correspond 
to the three key dimensionless parameters: |/3|, loq/uj and 
|a|. Region (1) is the zone where the formula for the in- 
dividual revival signals (|39[) derived within the adiabatic 
approximation is valid, and to summarize, we list here 
the restrictions on validity: 

• luq < 0.25a;: The adiabatic approximation is valid. 

• Ofi ^ /3 2 : Necessary for the eigenvalue and eigen- 
function simplifications leading to (|30[) . 

• |a| 3> 1: Validates the assumptions made in the 
Appendix concerning Sk(t). 

• 1/3 1 < 0.2: Imposed according to current experi- 
mental realizability. It is necessary for nearly de- 
generate states to have the same value of oscillator 
excitation number, N, and also for the simplifica- 
tions that lead to (150)1 . 

• | a/3 1 -C 1: Necessary for restricting the power series 
expansion of the Laguerre polynomial, to the 
first three terms (see (fS]3j). 




)/CO 



FIG. 8: Region (1): Parameter regime where the analytic 
formula for the collapse and revival dynamics derived within 
the adiabatic approximation is valid. Region (2): Parameter 
regime where the eigenspectrum derived within the adiabatic 
approximation is valid. Region (3): Parameter regime where 
the RWA is valid. 



For |/3| < 0.25, region (2) in Fig. [8] corresponds to the 
parameter regime where the eigen-spectrum of the sys- 
tem derived within the adiabatic approximation is valid. 
In region (3) of Fig. [5] we show the regime where the an- 
alytic formula for the collapse and revival si gna ls of two- 
qubit TC model derived within the RWA [35|] is valid. 
At resonance, the RWA results are valid for coupling 
strength as big as |/3| = 0.2. With the increase of de- 
tuning, the validity of the RWA is restricted to lower 
values of the coupling strength. We see that regions (1) 
and (3) are completely disjoint and thus the dynamics 
predicted within the adiabatic approximation cannot be 
derived from RWA calculations. 



IX. CONCLUSION 

In this report we extended the Tavis-Cummings model 
for multi-qubit interaction with a common oscillator be- 
yond its familiar RWA limits, in order to acomplish two 
goals: (a) to analyze two-qubit dynamics in the quasi- 
degenerate regime, following seminal work for a single 
qubit and oscillator by Irish [2l|, [22| others [H, [32j , 
and (b) to extend studies in this domain to include the 
dynamical behavior of quantum coherence in the form 
of qubit-qubit entanglement. We restricted attention to 
the regime where the qubit frequencies are much smaller 
than the oscillator frequency and the interaction coupling 
energy between the qubits and the oscillator is allowed to 
be an appreciable fraction of an oscillator energy quan- 
tum. 

We worked within the same adiabatic approximation 
introduced by Irish et al. [2l|. We showed that the RWA 
and quasi-degenerate regions do not overlap, and that 
the expressions derived lie completely outside the valid- 
ity regime of the RWA. We were able to compare fea- 
tures of single qubit dynamics with the corresponding 
extended results for two qubits, to identify features that 
are a consequence of having multiple qubits in the system 
and are qualitatively different from the single qubit case. 
An example occurs in the probability for two qubits to 
remain in their initial joint state. It is found to have a 
two-frequency beat note in the excitation-revival signals. 
This is absent in the single-qubit case. 

In cases of coherent-state preparation we were able to 
obtain convenient analytic formulas not previously avail- 
able, and showed that the analytic predictions compare 
favorably with exact numerical results. Expressions for 
individual collapse and revival signals were derived, pro- 
viding formulas for the width, height and time of indi- 
vidual revivals. 

Tracking of entanglement evolution, as a principal 
measure of intrinsically quantum coherence, is compli- 
cated even in the two-qubit case because several varieties 
of entanglement are present. We concentrated on qubit- 
qubit entanglement as our primary case study, which re- 
quired a trace over the oscillator's degrees of freedom. 
The 4x4 two-qubit density matrix yielded a compact 
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concurrence formula that included sequences of collapse 
and revival signals, and also indicated a route to control 
over them. Analytic expressions not previously available 
were derived for revival strength and timing. 

In this report, we have assumed that the system 
dynamics is not affected by interaction with a larger en- 
vironment. It is an interesting question to determine in 
what ways an external environment will severely or not 
severely impact the results presented. In order to lift this 
limitation, further analysis of the decoherence behavior 
of qubit-oscillator systems in the quasi-degenerate ultra- 
strong coupling regime is necessary. The possibility of 
generating non-classical states of the oscillator by letting 
it interact with multiple qubits in the quasi-degenerate 
ultra-strong regime is important but is not addressed in 
this report. 
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Appendix: CALCULATION OF S(t,uo) 
The infinite sum that we want to calculate is: 

oo 

S(t,u ) = P ( N ) cos (MNi\N )t), (A.l) 

N=0 

The Poisson distribution P(N) has an average and vari- 
ance equal to |a| 2 . If |or| is big enough, one can ap- 
proximate P(N) to a Gaussian with the same mean and 
variance and justify the replacement: 



P(N) 



\2N 



(A.2) 



An analytic form for the infinite sum is desirable but chal- 
lenging because of the Laguerre polynomial appearing in 
the definition of (Ni\N ). One notes that if \a/3\ <C 1, 
one can approximate the Laguerre polynomial by retain- 
ing only the first three terms of it to get: 

wo^lJVb) « cooe-*/ 2 (l - Nx + W^l gA (A.3) 

where x — (3 2 . When this approximation is justified we 
can insert it in summation (|A.1[) and get: 



S(t,uo) = ReY, P(N) 



N=0 



x exp 



IT 1 - Nx + 



N(N — 1) 



where we have defined a dimensionless scaled time by 
r = uj Q te- x/2 . 

Now we use the Poisson sum formula, according to which 
we get: 



S(t, ujq) = Re 

where 

S k (t,w ) = 



1 



w J 



P(0) exp (it) 



,k— — c 



(A.4) 



'0 

exp 



dnP{n)e 2i7rkn 
it I 1 — nx + 



n(n — 1) 



Using the replacement (IA.2I) . we see that S).(t, (Jo) be- 
comes a Gaussian integral, and when the excitation num- 
ber of the oscillator is great enough, \a\ 2 ^> 1, one can 
extend the lower limit of the integral to n = — oo and 
evaluate the integral analytically. The result is: 



S k {t,uj ) = 
where 



■exp($R e + i$r TO ), (A.5) 



Re 



i - (y + - 27rfc) 2 



4>; 



2(1 + (y//2)2) 

+ yf(y + yx/A-27rk))-\a\ 2 /2, (A.6) 

((l-(y + yx/4-2irk) 2 )yf/2 

, (A.7) 



2(l + (y//2)2, 

-2(y + yx/4- 2nk)) - 9/2 



and we defined: 



y = tx, 

f = \a\ 2 x, 

6 = tan- 1 (nkf) 2 . 



(A., 



The contribution of Sk(t,uio) to the sum S(t, lvq) will 
be maximum when <I>/j e is maximum, which occurs at 
times around y = 2itk. With this observation, and using 
i« 1, / <C 1, we can simplify the expression for <&R e by 
neglecting terms that are of the order of yx, (y — 2nk)f 2 , 
/ 3 , (y — 2nk) 2 f and higher powers of these terms to get: 



(y - 27rfc(l + f/2)) 2 . (A.9) 



2(l + (vrfc/) 2 ) 

Similarly, one can simplify the expression for $/ m to get: 

tan~ 1 (-Kkf) 2 1 , , , , 
*im = Y 12 - + ~ (y(l + /) - 27rfc/) . (A.10) 



Similarly, one can approximate the prefactor in (|A.5[) to 
be: 



1 



1 



l + (y//2)2 l + (7rfc/)2- 



(All) 
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Thus, Sk(t,LUo) takes the simplified form: 

s/ \ f -{T-T k ) 2 \a\ 2 x 2 >JR 

1 



(A.12) 



(l + (7rfc/) 2 ) 1/4 ' 
where <I>/ m is given by (|A.10I) and r k is defined as 

r k = 2nk(l + f/2)/x. (A.13) 



We now note that the amplitude of Sk(t,oj ) is much 
greater than P(0) = e - '"' . Thus, we can neglect the 
term |P(0) exp (it) from (|A.4I) . We also note from (|A.9|) . 
that for positive time r, we must have k > 0. Thus the 
expression for S(t,LUo) becomes: 



S(t, wq) = -Re 



oo 

E 

,fe=0 



S k (t,u>o) 



(A.14) 
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